library(rjags)
library(tidyverse)
library(ggmcmc)
library(kableExtra)
library(magrittr)
library(readxl)
if(!exists("i")){
  i <- commandArgs(trailingOnly=TRUE) %>%
    as.numeric
}
outdir <- file.path("..", "output", "timing.R")
if(!dir.exists(outdir)) dir.create(outdir)
set.seed(i)
mutations.per.year <- seq(1, 10, by=0.01)
mutations.per.year <- mutations.per.year[i]


unixFriendly <- function(x){
  x <- colnames(x)
  x <- tolower(x)
  x <- gsub(" ", "_", x)
  x
}
fname <- file.path("..", "data",
                   "Timing_Metrics.xlsx")
dat <- read_excel(fname,
                  sheet=2) %>%
  "["(-nrow(.), ) %>%
  set_colnames(unixFriendly(.)) %>%
  mutate(additional=mutations_ca-mutations_hg_ipmn) %>%
  filter(case != "MTP19")
dat
ipmn.age <- read_excel(fname,
                       sheet=1) %>%
  "["(1:2, ) %>%
  set_colnames(c("age", "n", "url")) %>%
  mutate(age=as.numeric(age),
         n=as.integer(n))
ipmn.age
ipmn.age2 <- ipmn.age[2, ]
pdac.age <- read_excel(fname,
                       sheet=1) %>%
  "["(5:6, ) %>%
  set_colnames(c("age", "n", "url")) %>%
  mutate(age=as.numeric(age),
         n=as.integer(n))
pdac.age
Tavg <- pdac.age$age[2] - ipmn.age2$age



y <- dat$additional
m1 <- jags.model("clock.jag",
                data=list(y=y, N = length(y), mu=mutations.per.year),
                inits=list(T=rep(5, length(y))),
                n.adapt = 2000, n.chains = 3)
jags_samples <- coda.samples(m1,
                             variable.names="T",
                             n.iter=5000,
                             thin=50)
jags_df <- ggs(jags_samples)
summaries<- jags_df %>%
  "["(grep("T", .$Parameter), ) %>%
  group_by(Parameter) %>%
  summarize(lower=quantile(value, 0.05),
            median=quantile(value, 0.5),
            mean=mean(value),
            upper=quantile(value, 0.95)) %>%
  mutate(mu=mutations.per.year)
saveRDS(summaries, file=file.path(outdir,
                                  paste0("timing_", i, ".rds")))


